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We consider the magnetic-field dependent spatial magnetization pattern around a general impu- 
rity embedded in a Heisenberg antiferromagnet using both an analytical and a numerical spin wave 
approach. The results are compared to quantum Monte Carlo simulations. The decay of the magne- 
tization pattern away from the impurity follows a universal form which reflects the properties of the 
pure antiferromagnetic Heisenberg model. Only the overall magnitude of the induced magnetization 
depends also on the size of the impurity spin and the impurity coupling. 
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I. INTRODUCTION 

The local magnetization around impurities in anti- 
ferromagnets have been studied by Nuclear Magnetic 
Resonance (NMR) experiments already since the early 
1970'sPEl The analysis of local Knight shifts has 
been expanded after the discovery of high tempera- 
ture superconductivity.^ Typically, the strongly corre- 
lated state is reflected by the observation of large al- 
ternating magnetic moments around static impurities^ 
which become especially strong in one-dimensionP An- 
other remarkable experimental tool is given by scanning 
tunneling microscopy (STM)P which offers the unique 
possibility of studying materials directly on the atomic 
scale. In particular, by coating the STM-tips with differ- 
ent magnetic materials,^ so called spin polarized scanning 
tunneling microscopy (SP-STM) has made it possible to 
study the magnetization of individual at orris P 

From the theoretical point of view, antiferromagnets 
are often represented by the isotropic Heisenberg model 
with static impurities. In this case the pinning of the 
order is a result of an interplay of the applied uniform 
magnetic field with impurities. The first theoretical stud- 
ies of impurities in an antiferromagnet date back to the 
1960's. 8 9 More recent research has made much progress 
in the und erstand ing of the impurity behavior in one- 
dimensiona l 4 * 10 * 1 ^ and two-dimensionaP^^' Heisenberg 
antiferromagnets. In particular, the magnetic response 
around a vacancy in an isotropic antiferromagnet was 
studied in Ref. |T5] using a hydrodynamic approach. In 
this work, we now extend those studies by considering 
the local magnetization using spin wave theory for a more 
general impurity type, which is given by a spin-So cou- 
pled to the host antiferromagnet with a general coupling 
Jo. One main result is that the decay constant of the 
magnetization is to leading order governed by properties 
of the host magnet, while the overall magnitude is gov- 
erned by properties of the impurity and its coupling to 
the host antiferromagnet. We complement our analytical 
spin-wave analysis with Quantum Monte Carlo (QMC) 
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FIG. 1: The canted spin state for classical spins. 6% £ [0, tt/2] 
is the angle between the spin i and a line drawn perpendic- 
ular to the applied magnetic field B. The angle fa £ [0, 2n\ 
parametrizes how much the spin i is rotated (a full rotation 
is indicated by the ellipse) about the applied magnetic field. 

simulations as well as a numerical spin wave approach for 
the case of calculating the magnetization on and close to 
the impurity site. 

II. HAMILTONIAN 

We consider the following Hamiltonian of a 
Heisenberg-type magnet in a magnetic field 

on a hyper cubic lattice where each site has Z nearest 
neighbors. We will start out with general site dependent 
couplings Jij and magnetic fields Bi and later specialize 
to the case of a single impurity in an otherwise uniform 
antiferromagnet in a homogeneous field. 

In order to treat the non-homogeneous Hamiltonian 
in Eq. (PD) with spin wave theory, let us first review in 
detail how to derive the expansion in fluctuations about 
an ordered classical state. The classical state of an an- 
tiferromagnet in a magnetic field is that of canted spins 
pointing partly along the z-axis, see Fig. [I] In order to 
parametrize this state we introduce rotated spins S' so 
that S[ z points along a direction parametrized by the 
angles 8i and , see Fig. [T] 
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The rotated spin components S' are related to the spin 
components in Eq. as 

Sf = (Sf sin 9i - Sf cos 6i) cos 4 - Sf sin fa 

Sf = (Sf sinOi-Sl* cosOJsmfc-S'Vcosfc (2) 

Sf = Sf cos 0, + Sf sinfli. 

Inserting these into Eq. ([!]) we get the Hamiltonian ex- 
pressed in terms of rotated spins for arbitrary angles, 
which will be determined later. In order to express the 
fluctuations about the ordered state we use the Holstein- 
Primakoff transformation^ on the rotated spins into 
bosonic operators 



S 



S'r = 




(3) 



Xcii - a]a]cLi ) + 



where expanding the square roots and using Sf = S' x ± 
iS' y yields 

ST = \[y ( a * + a « ~ i {^ a * a * + a ^ a ) + • ■ ■) (4) 

By inserting these expressions for 5" into the Hamilto- 
nian Eq. ([I]) we get terms H n with different powers n of 
bosonic operators. 

The zeroth order term in boson operators corresponds 
to the energy of classical spins oriented along the S' z 
axes. This is so because in the classical limit Sj — > oo 
the S' x and S' y components are overwhelmed by the S' z 
component which is proportional to S. The zeroth order 
terms read 

H Q = 2^ JijSiSj (cos 6i cos 9j cos((/)ij) + sin d t sin 9j) 



(ij) 



(5) 



where 



Because of the U(l) symmetry of 



spin rotations about the magnetic field axis Ho depends 
on the relative angles <pij . Minimizing with respect to <f>ij 
gives the condition 



JijSiSj cos( 



= or 7T. 



sin 







(6) 



For this to be a minimum 
-Jij cos(0y ) > 0, which means 
7r for an antiferromagnctic coupling and for a 



meaning that <f>ij 
of the energy one needs 
that d> 



ferromagnetic one. Equivalently — cos((j>ij) — Jij/\J%j\ = 
Vij. We will in the following select the rotation angle 
d>o so that it is either or w. With this choice, and the 
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FIG. 2: Couplings. Dashed lines indicate the coupling Jo to 
the impurity site (empty circle) while solid lines indicate J. 



minimization condition fcj — or 7r, all terms with sin < 
will be zero. Then the Hamiltonian can be written 



H 



cos 



) (s'fs: 



v tj S?S? 



<ij> 



-v^S 1 ? + sin(0 4 + Vii 9i) (vijS'fS'f + S' Z S?)} 
-J2Bi(S'* cosOi + S? siney. (7) 
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We will now specialize to the case of a single impurity 
embedded in an otherwise uniform antiferromagnet of 
spin-5 spins. We label the impurity site i — and allow 
for an impurity spin Sq which in general can be different 
from S. We take all bonds not connected to the impu- 
rity to be antiferromagnetic with a magnitude J. The 
bonds connected to the impurity are also equal, but of a 
different magnitude Jo and can be either ferromagnetic 
or antiferromagnetic, see Fig. [2j u denotes the sign of 
Jo- This antiferromagnet is placed in a magnetic field 
oriented along the z-direction with magnitude B. We 
have absorbed the Zeeman coupling into the magnitude 
of the magnetic field. In order to allow for a different 
gyromagnetic factor of the impurity spin and thus a dif- 
ferent Zeeman coupling, we label the magnitude of the 
effective magnetic field on the impurity site B which in 
general can be different from B. 

In order to simplify Eq. ^ we use an initial rotated 
frame that is given by a site independent value of 0i = 
9 for all sites i away from the impurity site to zeroth 
order. We will later allow for a site-dependent shift of 
9 in order to calculate the non-trivial local variation of 
the magnetization. For the impurity site i = we keep 
a separate angle 9 Q . Performing this ansatz the zeroth 
order term in boson operators takes the form 



H = -NS I cos 26» 



B sin 9 

+ZS {JScos29- \J \Socos(9 + v 9 )J (8) 
+BSsin9 - BoS sm9 . 

Minimizing this with respect to 9 and 9o in the thermo- 
dynamic limit, TV — > oo, determines the angles 9 and 9q 



sin 9 



B 



2SZJ 



(9) 
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and 



tan#o 



B 



JnlSZcOSt 



vq tan / 



(10) 



The zeroth order condition on 9 is identical to the one 
found for a uniform antiferromagnet in a homogeneous 
field and does not depend on the impurity. This is a 
natural consequence of taking a site-independent ansatz 
in the thermodynamic limit. 

When using the value of 6 obtained from Eq. ^ the 
terms that are of linear order in boson operators con- 
nected to the bulk behavior vanish. After also using the 
condition Eq. ( 10 1 only linear terms of bosons around the 



impurity are left 



H = 



-E 

<0j> 



(11) 



where the sum is restricted to run over the nearest neigh- 
bors of the impurity spin. This expression can be inter- 
preted as a local effective field in the rotated frame acting 
on the spins that are coupled to the impurity spin, which 
will cause a shift of the angles 9 over an extended range 
as we will see later. 

The constant C is given by 

C = J Q S Q zJ^u sin(u e + 9)- JSzJ^ sin29 (12) 
or equivalently when we use the minimization conditions 

° = vf (y^ B oCos0 o -Bcos^ . (13) 

The linear terms can also be written in terms of Fourier 
transforms 



at = — ;= \ are 

k 



ik-ri 



as 



C V / 



(14) 



(15) 



known from standard spin- wave theory.^ The neglected 
quadratic impurity terms can in principle lead to a renor- 
malization of the overall magnitude in the local order 
around the impurity. However, this effect is known to 
be surprisingly small from numerical studies,^ so that 
we can omit those terms for now in order to calculate 
the magnetization around the impurity. We will include 
them later when considering the magnetization of the 
impurity spin itself. 

The quadratic term can be diagonalized by the canon- 
ical transformation 



which results in the quadratic Hamiltonian 

2^ = £"^s + ?£te-^) 



where 



Al- 
ls 



JSZ 



Bl which becomes 

k 

^(1-7*:) (l + co S 20 7jE ) 



(17) 



(18) 



(19) 



The transformation coefficients obey ul — vi = 1, ul + 

v l = ^s/ u S and 2u k v k = - B k/"k- 

Using the quadratic bulk Hamiltonian we can calculate 
the following expectation values 



A = (ajOj 



m = {a 



k 

J\J ^ / ^k U k V k 
k 



(20) 



for nearest neighbor sites i and j. Note that the bulk 
nature of the quadratic term dictates that these expres- 
sions do not depend on i and j. At this stage we trun- 
cate higher order terms in the Hamiltonian. Therefore we 
have reduced the problem to a solvable bulk Hamiltonian 
in Eq. (161 together with an impurity term in Eq. (15). 



where we have defined jk = 2(cosfc K + cosfc^ + . ..)/Z 
where the k's are given in units of the inverse lattice 
spacing. 

For the quadratic terms we will as a first approxima- 
tion keep only the terms that are leading order in N. 
Therefore, the quadratic terms are identical to those in 
the absence of an impurity 

^2 Ulk = 2 E { A AH + B i a ^-k + h-o) (16) 



where A R = JSZ{cos 29 - 7 £ sin 2 6) + B sin 9 = JSZ(1 - 
7£ sin 2 9) and Br — JSZ cos 2 6*7g which are also 



III. MAGNETIZATION AWAY FROM THE 
IMPURITY 

The magnetization in the direction of the field Mf = 

(SI) is 

M? = (Sf ) cos Bi + (Sf ) sin 0i. (21) 

Expressed in terms of bosons the above expression is up 
to quadratic order 



Ml RJ sin I 



cos 




(22) 
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To calculate these expectation values in the presence of 
the impurity we perform a shift of the boson operators 

di -> a, + on (23) 

so as to get rid of the remaining linear terms in the Hamil- 
tonian in Eq. ( 15 ). This is equivalent to a site dependent 
variation of the angle 9i. The impurity induced shift is 
given by 



N^A % + B % 

k 



For future convenience we parametrize 



(24) 



(25) 



in terms of constants / and g which to leading order in 
1/S are obtained from Eq. ( 16 1; / = JSZ and g — cos 29. 



Shifting the boson operators gives the following expres- 
sion for the magnetization 



Mf ss sin 9, 



(Si 



■ cos 9i (a* + on) . 

(26) 



Since the shift of the boson operators has eliminated the 
linear terms, we can now use the usual bulk theory to 
calculate the corresponding expectation value n 
in Eq. (20 1. Thus the magnetization takes the form 



{aUi) 



Mf k, sin (9 (S- \ ai \ 2 - n) + \j - cos 9 {a* + on) ,i ^ 0. 

As is shown in the Appendix, on is real and changes 
sign depending on which sublattice i belongs to with 
e l & p = (-l)*<+v«+2 4 where Q = (tt,tt,...) is the anti- 
ferromagnetic wave vector. Therefore, it is convenient 
to write on — (— l) Xi+Vi+Zi cti and to divide the magne- 
tization into an alternating and a non-alternating part. 
Using the assumption that on does not vary rapidly, the 
alternating(non-alternating) magnetization on site i is 
obtained by taking half of the magnetization on an odd 
sublattice site i and subtract (add) half of the magnetiza- 
tion on the neighboring even sublattice sites surrounding 
site i. Therefore, the non- alternating part takes the form 



nalt,2 



= sin i 



(S- ri- 



al) 



(28) 



which will decay rapidly to its uniform bulk value. This 
non-alternating part is not our primary focus here. In- 
stead we will focus on the alternating part which does 
not decay as rapidly. To leading order the alternating 
magnetization is 



M. 



alt, i 



= — v25cos# OLi 



(29) 



thus b\i dictates its behavior. The sum in Eq. ( 24 1 can be 



carried out by expanding the integrand about the min- 
imum of the denominator which is at the antiferromag- 
netic point Q = (jr, tt, . . .) as shown in the Appendix. 
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FIG. 3: (color online) M a z „ vs. distance from the impurity r 
on the square lattice. The circles are quantum Monte Carlo 
data while the dashed line (red) is a plot of the analytic result 
Eq. (291 using g = cos 2(9. The result where we have taken 



into account 1 / S'-corrections for + is shown as the solid 
line (green). Here S = S = 1/2, Z = 4, B = B = 0.4 J and 
Jo = 0.1J. 



Carrying out this expansion for the case i ^ 0, we get in 
D = 2 and D = 3 dimensions 



CZ 



Mn/d), 
e- r */ d /(2n), 



D 



,i^0. 



(30) 



D = 3 



where Ti — \J x\ + yf + zf is the distance from the impu- 
rity in units of the lattice spacing and Kq is the zeroth 
order modified Bessel function of the second kind which 
decays as e~ ri ^ d /y/ri for large arguments. The charac- 
teristic decay scale is 



d = 



<J 



z{i-g) 



(31) 



in both cases. The result in Eq. (30) is the main result of 



this section for the induced magnetization by the general 
impurity model, which will be compared to Monte Carlo 
results in the following. Note, that the shape and the de- 
cay scale d is universal and only depends on properties of 
the host magnet in the bulk. Only the constant prefactor 
C in Eq. (13) depends on impurity properties S^, an d 
Bq. With the expression g = cos 29, the decay constant 
is d= [cos2<V(2Zsin 2 (9)] 1 / 2 . 

In Fig. [3] we have plotted a comparison of M* lt cal 



culated using the expression in Eqs. ( 29 )-( 30 ) and re- 
sults from a QMC simulation. The QMC simulations 
were carried out using the stochastic series expansion 
technique^ using directed-loop updates^ at a low tem- 
perature T/J — 0.05 on a 128 x 128 square lattice. As 
can be seen from Fig. [3] the leading order analytical result 
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FIG. 4: (color online) C/J vs. impurity coupling Jo for im- 
purity spin So = 1/2 (upper panel) and So = 1 (lower panel) 
for different values of the magnetic field B/J indicated by the 
numbers above each curve on the left side. Here S =1/2 and 
Z = 4, B = B. 



decays faster than the QMC result. However the decay 
d depends crucially on the exact expression for + B^ 
which we have approximated with its leading order value 
d = [cos26»/(2Zsin 2 6>)] 1/2 . In fact, we can do better by 
including l/S corrections. Taking into account l/S cor- 
rections to A^ + B£ and to the angle sin (9, we get 



A k+ B k 



JSZ 



1 



sin 2 6»- 



+7fe cos 29 



2n + 2A + m 

2s ' "" 2.s 

2n + 2m + 2A + S 



2s 



sin 2 (9- 



,2n + 2m + 2A-5 



(32) 



This result can also be inferred from Ref. 17, The l/S 
corrections give modified expressions for the constants 
/ and g, which leads to a better agreement with the 
Monte Carlo data in Fig. [3] By allowing also another 
classical angle 9\ for the impurity nearest neighbor spins 
the agreement with QMC close to the impurity site can 
be improved at the expense of having more complicated 
analytic expressions. To connect our result in Eqs. ( 29 1- 



( 30 1 to that obtained in Ref. [TS] for the induced mag- 
netization around a vacancy (Jo = 0) we observe that 
for k close to Q but \k - Q\ > [8sin 2 6»/cos26»] 1 / 2 the 
dispersion Eq. ( 19 1 is linear with a spin- wave velocity 
c = 2 JSV2 cos 20. In the limit B — > this becomes the 
well-known leading order spin wave theory result for the 
spin wave velocity of an antiferromagnet. Combining this 
with Eq. ([9]) we see that the decay constant of Ref. [TBI be- 
comes c/B = [cos 28/(8 sin 2 0)] 1 / 2 , which equals the lead- 



FIG. 5: (color online) C/J vs. magnetic field B/J for different 
values of the impurity spin and coupling denoted by (So, Jo). 
Here B = B, S = 1/2 and Z = 4. 



ing order result for the decay constant d. Similarly, we 
can compare the factor multiplying the Bessel-function 
Kq. In the case of a vacancy Jo = our expression for 
C = —(S/2) 1 l 2 BcosO so that the prefactor becomes 



'2SCOS0- 



C 



^fg 2 



B 
2^~J 



(33) 



where we have used / = JSZ and g = cos 26 and approx- 
imated cos 9 f=a 1 which is valid for low magnetic fields. 
This is to be compared to the expression m max SB/(2'?rp s ) 
obtained in Ref. 15, When inserting the leading order ex- 
pression m max = 5, p s — JS 2 we see that the two results 
become equal. 

For larger fields the use of the renormalized zero field 
spin- wave velocity c in Ref ll5l is not so natural, however. 
As the decay depends heavily on the behavior of A^ + B? 

around k — Q where the dispersion is quadratic in a finite 
field, it is more natural to relate the decay constant to 
the effective mass of this minimum. For finite but not 
too large fields the dispersion around Q can be written 
Jt- where the effective mass is ™ — 2Zsin 9 



B cos 29 ■ 

It is then straightforward to see that the l eadin g order 
decay constant also can be written d = 1/vBm. 

While the decay of the induced alternating magnetiza- 
tion pattern is governed by the properties of the uniform 
magnet, the magnitude of the alternating magnetization 
is given in terms of the prefactor C in Eq. ( |13[ ) , which de- 
pends on impurity properties as shown in Figs. [4] and [5] 
For impurity spin So — 1/2 and coupling < Jo < 1 , the 
prefactor C is negative and rather small. For Jo = J it 
vanishes completely because it corresponds to the uni- 
form case. For ferromagnetic couplings Jo < 0, \C\ 
gets larger with increasing magnetic field B/J. Thus we 
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FIG. 6: Orientations of the magnetization close to the im- 
purity. The impurity spin is the middle circle. Open circles 
indicate that the magnetization is pointing along the applied 
magnetic field while filled circles indicate the opposite orien- 
tation, a) C < 0. b) C > 0. 



expect a substantial induced alternating magnetization 
pattern for ferromagnetically coupled impurities. Note, 
however, that when the field gets larger the magnetiza- 
tion pattern decays faster with distance from the impu- 
rity. For an So = 1 impurity, \C\ is no longer necessarily 
small for antifcrromagnetic couplings and it changes sign 
at a small positive value of Jo/ J. The sign change sig- 
nals a sublattice change in the magnetization pattern as 
indicated in Fig. [6j where for a ferromagnetic impurity 
the magnetization follows the pattern shown in Fig. [6] a). 
This pattern extends also to weak antiferromagnetic cou- 
plings up to a critical value of Jo that depends on the 
magnetic field where it becomes favorable to interchange 
the orientation of magnetization on the two sublattices 
while keeping the impurity spin oriented along the field. 
This results in the pattern shown in Fig.[6]b). For large 
values of Bj J and for all couplings except large antiferro- 
magnetic ones, \C\ increases linearly with field strength 
Bj J as shown in Fig. [5] For Sq = 1 and a small antifer- 
romagnetic coupling Jo, C changes sign as the magnetic 
field is increased, second curve from the top in Fig. [5] 
Thus a change in the sublattice rearrangement in Fig. [6] 
can also happen for a fixed Jo as the magnetic field is 
varied. The exact point where C reverses sign is special, 
because when C = the spin-1 impurity appears to have 
no effect on the host spins of the surrounding antiferro- 
magnct. Therefore, the field and/or the coupling can be 
tuned in such a way that the impurity becomes almost 
invisible to the bulk, i.e. very little scattering occurs. 



IV. MAGNETIZATION OF THE IMPURITY 
SPIN 

At the impurity site the leading order magnetization 
is obtained by the classical expression 



Mn = Sb sin I 



'o- 



(34) 



For So = 1/2 and J > this gives a reasonable agree- 
ment with the QMC data, as is seen in Fig. [7] However 
for other spins and ferromagnetic couplings Jo < the 
result is rather far of the QMC result. Thus it is neces- 
sary to also take into account the quantum corrections 




FIG. 7: (color online) Magnetization at the impurity site for a 
spin- 1/2 impurity coupled to a bulk spin- 1/2 antiferromagnet 
by a coupling Jo. The filled black circles are results from 
Quantum Monte Carlo simulations. The dashed line (green) 
is the classical result coming from Eq. (341, and the solid line 
(red) is the numerical spin wave result. 



to Eq. (34). However, these quantum corrections are dif- 
ficult to calculate analytically. This is because for the 
impurity itself it is necessary to include explicitly the bi- 
linear terms connecting the impurity site to its neighbors 
in addition to the quadratic bulk part in Eq. ( 18 ). These 



impurity terms induce non-local interactions in fc-space, 
thus an analytic diagonalization becomes difficult. In 
order to solve this we will instead numerically diagonal- 
ize the quadratic boson Hamiltonian as described below, 
which gives much better results shown in Fig. [7| As this 
method is numerical there is no need for the restriction 
of keeping only two angles 9q and 9. Thus we will instead 
keep track of all the angles 9i. This has the consequence 
that all linear boson terms vanish when using the values 
of the angles obtained from minimizing the zeroth order 
term, as will be shown below. 

As a function of all angles 9i the zeroth order term is 

H Q = ^ -\ J ij\ S i S j C0S(6> 4 + Vij9j) - ^2 B i S i Sm 

<ij> i 

(35) 

where we have used the minimization condition for the 
</>'s. Minimizing Hq with respect to 9i we find 

^2 I J *:> I Sj sin {9, + v tj Bj) - Bi cos 0< = (36) 

where the sum is restricted to run over the nearest neigh- 
bors ei of site i. This condition is equivalent to the equa- 
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tion 



tan 6*, = 



Ej= e( \Jij\Sj COS 6 j 



(37) 



The operators S'fS'f, S'fS'f and the magnetic field term 
in Eq. ^ give the linear terms of the Hamiltonian 
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FIG. 8: Geometry of a 6 x 6-lattice. The open circles mark 
sites where the boundary condition is imposed. The filled 
circles are sites where the angles are being calculated. The 
small circle is the impurity site. Periodic boundary conditions 
are used. 



\Jij\Sj sin(0i + VijBj) - Bi cos^ . (38) 



By comparing this to Eq. ( 36 1 we see that the minimiza- 



tion of the constant terms leads to the vanishing of the 
linear terms. 

The quadratic terms are 



Ho 



E ( C0S (^ + V ij0j) ~ V ij) { a \ a 3 + a ] a i) 

<ij> 

+JijVij cos(6>i + VijOj) (^SjajcLi + Sia^a^j 
+Jij\J ^^{cosidi + v^O^ + Vij) (aiOtj + a| at) 



sm f,a„-a,- 



(39) 



which can be written in the form 

H 2 = ^2 ( a i^ij a j + Oi-AijOj + a \ B ij a ) + "-/',',"; ) + G 

(40) 



where the constants are 



3=&i 



(41) 




Aij = Jij v ' (cos(6»i + VijOj) - v^) 6 <i j > (42) 
— sin^ + ^2 '^VikSk cos(9i + v ik 9j) J 8 i3 - 



and 



\J Si Sj 



Bij = v " ' ~ J i '•()>; W, + VijOj) + v^) 5 <ij> (43) 

where 5 < ij > is 1 when i and j are nearest neighbors and 
zero otherwise. 



mg 



numerical 
This 



In order to numerically diagonalize Eq. (40) we will 
first find the 

Eq. m 



values of the 9^s by solv- 
is achieved by the relaxation 



method where the boundary condition is specified as 
sin ^boundary = B/2SZJ and an initial guess for the an- 
gles on other sites is made as indicated in Fig. [8j Then 
the lattice is traversed site by site and new angles are 
computed using Eq. (37). This step is repeated until 



convergence. It is known that this procedure converges 
slowly. However for typical lattice sizes (28 x 28) used 
here this is not an issue of practical importance. Having 
determined the angles numerically we proceed to diago- 
nalize the quadratic Hamiltonian. 

We begin by forming the 2N column vector a = 
(ai, a2, . . . , a,N, a\, 0,2, ■ ■ ■ , ajv) T where we have numbered 
the lattice sites in a consecutive fashion from 1 through 
N. The components of a obey the following commutation 



relation 



CU, 



= f]ij where r) = ( ljv x 



o 

-ljv X .V 



With 



this notation the quadratic Hamiltonian takes the form 



H = a f 2)o 



(44) 



where 2) is the 2N x 2iV-matrix with entries from the 
quadratic part of the Hamiltonian 



D 



A B 

B* A* 



(45) 



We seek a 2N x 27V Bogoliubov transformation matrix t 
that transforms a into new bosonic operators b: a = tb. 
In order for the entries of b to obey bosonic commutation 
rules the matrix t must obey 



(46) 



Inserting o = tb into the Hamiltonian (44) we seek a t 



that fulfills the commutation condition Eq. ( 46 ) and that 
makes tt Di = € where £ is diagonal. However it is not 
always possible to find such a diagonal matrix. When 
the Hamiltonian contains zero modes associated with a 
continuous spectrum one will never be able to write the 
free particle operator p 2 as a b^b term alone. However 
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such a term can always be written as tfb + btf — bb — 
6 t 6 t with the proper rescaling of operators. Thus we will 
seek a matrix € that is almost diagonal in the sense that 
for massive modes it has only entries along the diagonal 
while the continuous parts of the spectrum is represented 
by Is or -Is in appropriate places. More specifically we 
are seeking a matrix t that makes t^2)t into a 2N x 2N- 
matrix (£ of the form 



£ = 



Z Z 

h J 

Oz o 2 
V Jg ij 



(47) 



where E e is a diagonal e x e matrix of positive energies 
which represents the discrete harmonic oscillator energies 
associated with e gapped modes. Here Z is a z x z-matrix 
of zeros that represents z proper zero modes where the 
harmonic oscillator energy is zero, I z and J z are describ- 
ing the z improper zero modes associated with a contin- 
uous free-particle spectrum, I z is a z x z-diagonal unit 
matrix, and J z is a z x z diagonal matrix with diagonal 
entries either +1 or — 1. The sign distinguishes between 
operators of the type x 2 and p 2 . Empty entries indicate 
zeros. The procedure of finding such a t is outlined in de- 
tails in Ref . [2T| We have implemented this on a computer 
and find that the procedure works very well. 

In the absence of linear terms the magnetization is 
given to quadratic order by 



(5?> =Bin0 4 - <atoi>) 



(48) 



The value of sin 6?, is known from the minimization of 
the classical term, and (ajaj) can be obtained from the 
transformation matrix t. Without loss of generality the 
matrix t can be written 



t = 



U V* 

V u* 



(49) 



where U and V are N x TV-matrices. Then the expecta- 
tion value {a\aj) is 



aUi) 



+u: J v* k {b]b{} + v lJ u tk {b j b k : 



(50) 



We will compute the expectation value in a state with 
low energy. For massive modes we pick the ground 
state to be the vacuum state and then only the sec- 
ond term contributes (bjbV) = Sjk- The situation is not 
so simple for the improper zero modes. An improper 
zero mode tfb + 66^ ± bb ± tftf can be written as the 
momentum squared operator 2p 2 (the minus sign) or 
the position squared operator 2x 2 (the plus sign) using 
b = (x + ip) and = 4= (x — ip) . Thus it is clear 
that its spectrum is continuous. 



For each improper zero mode we choose instead to com- 
pute the expectation value in a Gaussian stated charac- 
terized by a width w. Specifically 



tp(x) 



1 



1/4 



-\(x/w) 2 



(51) 



In this state the expectations values of the energies are 

1 



(p > = 2 W 
(x 2 ) = \w 2 



(52) 
(53) 



while the expectation values of the operators needed in 
(a\ai) are 

(6*6) = (w 2 +w-' 2 ~2) /4 (54) 
(66 f ) = (w 2 + w-' 2 + 2) /4 (55) 
(fttftt) = (bb) = (w 2 - w- 2 ) /4 (56) 

Using this the expectation value (aj a i) takes the form 



+ ^\U: ] -V lJ \ 2 -2(\U l3 \ 2 -\V tJ \ 2 )j(57) 

We will refer to the last sum in the above as the zero 
mode(s) contribution, and we have allowed for a separate 
width Wj for each improper zero mode. We will choose 
values of Wj so that the total energy of the improper zero 
modes is equal to that of the lowest finite energy mode. 
This choice is made to avoid divergences and at the same 
time still justify calling them zero energy modes. In our 
case, in the presence of a magnetic field, there is only one 
improper zero mode, and it turns out that the precise 
value of the w is not important quantitatively for the z- 
axis magnetization. In all cases we have looked at here, 
the zero mode contribution is negligible and we might as 
well neglect it completely. This is in contrast to the one 
dimensional case where the zero modes dominate and are 
responsible for the divergences of spin- wave theory in the 
infinite volume limit. 

The results from this numerical diagonalization on a 
28 x 28-lattice is shown in Fig. [7] for Sq = 1/2 alongside 
the classical result and results from QMC simulations for 
the square lattice at a fixed value of the magnetic field 
B/J = 0.4. Fig. [9] is similar but for = 1. One can 
see that the numerical diagonalization procedure com- 
pares much more favorably to the QMC data than the 
classical result does. Especially for antiferromagnetic Jo 
the agreement is very good. For large ferromagnetic Jo 
the agreement is worse which we believe is related to 
the truncation of the Hamiltonian at quadratic order in 
boson operators. The main feature of the curves is a 
maximum at Jo which reflects the trivial fact that an un- 
coupled (isolated) impurity will point along the magnetic 
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FIG. 9: (color online) Magnetization at the impurity site for 
a spin-1 impurity coupled to a bulk spin-1/2 antiferromagnet 
by a coupling Jo. The filled black circles are results from 
Quantum Monte Carlo simulations. The dashed line (green) 
is the classical result coming from Eq. (341, and the solid line 
(red) is the numerical spin wave result. 



field. In fact the impurity spin will point along the field 
for most couplings except very large antiferromagnetic Jo 
for So = 1/2. 

For sites in the neighborhood of the impurity we can 
also compare the analytic and the numerical spin wave 
calculation to the QMC results. In Fig. [lO]we show the 
magnetization for an So = 1 /2 impurity at different posi- 
tions (xi, yi = 0) close to the impurity. The different lines 
are for the various values of the impurity coupling Jo and 
the different symbols indicate the method used. In com- 
paring the methods we see that the analytic result lies 
reasonably close to the QMC data except for the nearest 
neighbor point where the numerical spin wave calcula- 
tion give a better approximation to the QMC data. For 
a fixed value of Jo one can see that the magnetization ex- 
hibits a predominantly alternating pattern with a mag- 
nitude that is largest for ferromagnetic couplings Jo < 
as predicted in Fig. [4] As the ferromagnetic coupling Jo 
becomes smaller the magnetization of the impurity spin 
increases, Fig. [7j while the surrounding pattern is not 
much affected. On the antiferromagnetic side, Jo > 0, 
the magnetization of the impurity spin decreases accom- 
panied also by a decrease in the amplitude of the mag- 
netization oscillation away from the impurity. At Jo = J 
the oscillation pattern vanishes completely. For strong 
antiferromagnetic couplings Jo > J there is almost no 
induced magnetization on the sites surrounding the im- 
purity, but the magnetization of the impurity spin be- 
comes smaller than the average magnetization and can 
even become negative for strong enough Jq. 



FIG. 10: (color online) Magnetization as a function of hor- 
izontal distance x, from the impurity site as calculated by 
QMC (solid circles), numerical spin waves (triangles) and the 
anaytic spin wave theory (squares). So = 1/2, S = 1/2 and 
B = B = 0.4 J. The colors are for different values of Jo/ J =; 
— 2(solid, black), 0(long dashed, red), 0.1 (dot-dashed, green) 
and 0.5 (dot-dashed, blue). QMC error bars are smaller than 
the size of the solid circles, and both the QMC and the nu- 
merical spin wave calculations are carried out on a 28 x 28- 
lattice. 



For the So = 1 impurity the magnetization pattern 
around the impurity is shown in Fig. |11| Again the os- 
cillations are large for ferromagnetic Jo- As Jo — > the 
magnetization of the impurity spin increases while the os- 
cillating pattern around it decreases. Then as Jo becomes 
antiferromagnetic the magnetization oscillations increase 
again, but now the sublattice pattern has changed to 
the pattern in Fig. [p]b), consistent with the fact that C 
changes sign in Fig!~|4] The amplitude of the oscillations 
saturates as Jq becomes even stronger. 



V. DISCUSSION 

We have presented results for the magnetization 
around a general impurity in a Heisenberg spin-S antifer- 
romagnet in a magnetic field. Away from the impurity 
we find that the induced magnetization is dominantly a 
staggered magnetization in the held direction. We have 
calculated this alternating magnetization, and our results 
are in reasonable agreement with extensive QMC simu- 
lations that we have also carried out. One important 
feature of the spin wave result is that the parameters of 
the impurity model only affect the overall prefactor C of 
the magnetization while the scale and shape of the decay 
are universal and only reflect the properties of the host 
magnet and the applied held. We have analyzed how the 
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FIG. 11: (color online) Magnetization as a function of hor- 
izontal distance Xi from the impurity site as calculated by 
QMC (solid circles), numerical spin waves (triangles) and the 
anaytic spin wave theory (squares). So = 1, S = 1/2 and 
B — Bo = 0.4 J. The colors are for different values of Jo/ J = 
— l(solid, black), 0(long dashed, red), 0.2 (dashed, green) and 
1 (dot-dashed, blue). QMC error bars are smaller than the 
size of the solid circles, and both the QMC and the numerical 
spin wave calculations are carried out on a 28 x 28- lattice. 



prefactor C depends on impurity properties and found 
that the effect on the alternating magnetization is largest 
for ferromagnetically coupled impurities and generally in- 
creases with magnetic field. In order to calculate the 
magnetization at the impurity site we have described in 
detail how to diagonalize the quadratic spin wave Hamil- 
tonian numerically. This approach agrees well with the 
QMC calculations and we have outlined how the mag- 
netization of the impurity spin depends on the coupling 
strength of the impurity to its neighbors. 

In summary the results can be used to predict the de- 
tailed local magnetization pattern around general mag- 
netic and non-magnetic impurities in isotropic antiferro- 
magnets, e.g. from doping Zn, Co and Ni in copper-oxide 
antiferromagnets. In most real materials the effects from 
crystal fields and other anisotropies are also important, 
but our calculations provide the first step, before other 
possible terms in the Hamiltonian are taken into account. 
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Appendix: Sum 



The sum 



N ^ ' 

k 



91k 



(A.l) 



for f can be written 



/ = 



1 1 + 91% ~ 1 
g N 4- 1 + ff7fc 

k 



This sum can be calculated by expanding the denomi- 
nator about the antiferromagnetic point Q = (7r, tt, . . .). 
Shifting the fc summation k — > k + Q and expanding the 
denominator to order fc 2 we get 



Ak-r 



9 N ^ 1- g + gkl/Z 



(A.3) 



where Z is the coordination number of the lattice. This 
can also be written 



N 



d 2 k 2 



(A.4) 



where d 



z {i-g) • ^ ne sum i s calculated by trans- 
forming it into an integral and using polar coordinates 



Ak-r 



-y 

N ^? 1 + d 2 k 2 27TC? 2 r/d 



K (r/d), D - 2 
, , (A.5) 

7(2r), D = 3 

where Kq is the zeroth order modified Bessel function of 
the second kind. Putting this together we get 



Ze iQ - 



K (r/d), D = 2 
e -r/d/(2r), L> = 3, 



(A.6) 



where e lQ ' ? = (-l) x *+Vi+^ 
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